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Numerical modeling of the diffusional transport associated with high-tempera- 
ture corrosion processes is reviewed. These corrosion processes include external 
scale formation and internal subscale formation during oxidation, coating deg- 
radation by oxidation and substrate interdiffusion, carburization, sulfidation 
and nitridation. The studies that are reviewed cover such complexities as con- 
centration-dependent diffusivities, cross-term effects in ternary alloys, and 
internal precipitation where several compounds of the same element may form 
(e.g., carbides of Cr ) or several compounds exist simultaneously (e.g., carbides 
containing varying amounts of Ni, Cr, Fe or Mo). In addition, the studies 
involve a variety of boundary conditions that vary with time and temperature . 
Finite-difference (F-D) techniques have been applied almost exclusively to 
model either the solute or corrodant transport in each of these studies. Hence r 
the paper first reviews the use of F-D techniques to develop solutions to the 
diffusion equations with various boundary conditions appropriate to high- 
temperature corrosion processes. The bulk of the paper then reviews various 
F-D modeling studies of diffusional transport associated with high-temperature 
corrosion. 


KEY WORDS: modeling; numerical modeling; numerical techniques; finite-difference tech- 
niques; oxidation; corrosion; carburization; nitridation; diffusion. 


INTRODUCTION 

In this paper, numerical modeling will be limited to the methods and tech- 
niques describing the diffusional transport within an oxide scale, subscale, 
or alloy associated with high-temperature corrosion processes. The diffusing 
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specie(s) may be moving inward as a corrodant. such as occurs during 
internal oxidation, sulfidation or carburization, or outward as a solute, such 
as occurs during selective oxidation to form an external protective scale. 
These diffusional processes are described by the well-known Pick’s laws of 
diffusion, given as 


and 



(First Law-) 


( 1 ) 


vt dX\ dXj 


(Second Law) 


(2a) 


where J is the flux and C is the concentration of the diffusing specie, D is 
the appropriate diffusion coefficient, and X and t refer to distance and time, 
respectively. In an attempt to maintain clarity, C and N will refer to the 
elements carbon and nitrogen, while C and N will refer to a concentration 
parameter and to the number of nodes in a grid, respectively. Often, D is 
assumed independent of concentration (D^j(C)) such that Eq. (2a) 
becomes 


dC ^ / d 2 C\ (Second Law) 

vt ~ \dX 2 ) (£>#( O) 


When D cannot be assumed independent of concentration. Eq. (2a) can be 
rewritten as 


dC d 2 C dD dC 

— D ■ + 

vt dx~ dx rx 

and by using the chain rule for dD/ dX, Eq. (2a) finally becomes 

( C d 2 C dD ( dC\ (Second Law ) 

= D — — (2c) 

dt cX 2 dCXdxJ (D=/(C)) 

In most high-temperature corrosion problems, one knows the tempera- 
ture (even if it is not constant) and desires to describe or predict the extent 
of the attack at a given time. Typically, a knowledge of the concentration 
profile (C vs. X ) for the specie(s) of interest (even if the specie is precipitated 
as an oxide, carbide, or sulfide), or the thickness of an internal subscale, is 
sufficient to describe the extent of attack. Throughout this paper, subscale 
will refer to a zone within an alloy containing a dispersion of one or more 
precipitates of oxides, carbides, sulfides, etc. Elence, a solution to Fick’s 
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laws, yielding concentrations and interface positions, is generally sufficient 
to yield the desired information about the extent of the corrosive attack. 

Analytical solutions to Fick’s laws are available to describe many situ- 
ations involving high-temperature corrosion. Examples of these solutions 
cover such areas as solute oxidation in binary alloys, 1 internal oxidation, 21 
moving boundary problems, 4 and solute transport during oxidation of ter- 
nary alloys. 5 In addition, generalized solutions have been compiled in stan- 
dardized texts such as those by Crank 6 and Jost. 7 However, each of these 
solutions requires some form of simplifying assumptions, such as constant 
diffusion coefficients, isothermal conditions, semi-infinite geometries, or 
time-independent boundary conditions. Assumptions such as these are often 
unacceptable for many real-world problems. 

The utility and usefulness of numerical solutions are their flexibility 
to accommodate these real-world complexities that occur in many high- 
temperature corrosion problems. These complexities include strongly 
concentration-dependent diffusion coefficients (e.g., D in f3 NiAl), time- 
varying boundary conditions (e.g., the discontinuous rate of solute consump- 
tion due to oxide spallation during cyclic oxidation), complex initial condi- 
tions (e.g., a near-surface solute enrichment due to a coating treatment or a 
solute depletion due to a preoxidation exposure), finite specimen geometries, 
precipitation of multiple compounds of one or more solute elements (e.g., 
carbide formation), and nonisothermal conditions affecting the diffusion 
coefficients, boundary conditions or solubility products. 

Various numerical techniques have been employed to solve the differen- 
tial equations associated with diffusional transport, however, the most com- 
mon technique has been that of finite differences (F-D). To the author’s 
knowledge, the F-D technique is the only numerical method which has been 
applied to modeling transport associated with high-temperature corrosion, 
although other techniques 8 could certainly be applied. Hence, only the F-D 
technique will be discussed in this paper. In the following sections, this 
technique will be described followed by various examples in the literature 
of how r F-D techniques have been used in the area of high-temperature 
corrosion modelling. 

BASICS OF THE F-D TECHNIQUE 

The initial step in the F-D technique is to establish a grid across a 
region of the material over which diffusion will occur. The grid consists of 
N nodes with a spacing between nodes of AX. Typically, the node spacing 
is uniform although a finer or coarser, grid spacing can be used in regions 
w here the concentration gradient is anticipated to be either steep, or shallow, 
respectively. Each node of the grid has a specific concentration associated 
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with it. Fick's laws (Eqs. (1) and (2)) are replaced with F-D equivalents 
based on small differences in concentration, AC, distance, AX, and time, A/, 
as shown in the following section. A solution yielding the concentration 
profile at some time t is derived by solving the appropriate F-D equivalents 
for small time increments (Af) in an iterative manner. These iterations, or 
time steps, are continued until the At increments sum to the desired time t. 
Because of the typically large number of repetitive and tedious calculations 
made each iteration, F-D solutions are ideally handled by a computer. 


Explicit and Implicit F-D Methods 

The F-D equivalents to Fick's laws are based on Taylor series 
expansions. 9 The first order, central-difference equation to Fick's first law 
(Eq. (1 )) is given as 


+ ' ' ) (n = 2toN-\) (3) 

V 2A3f / 

where the subscript refers to the node number and the superscript / refers 
to the time increment or iteration. In this section dealing with the F-D 
technique, the subscript to C ( 1, 2, 3, . . . , N) will refer to the node number, 
whereas in the following application section of the paper, the subscript to 
C (O, C, N, etc.) will refer to a corrodant or solute. The F-D equivalent of 
Fick’s second law' can be given by either an explicit or implicit representation. 
The explicit form for the case where D^f(C)) is given as 

c:\-aj /CH.-2CK’;, 

At ~ 2 \ (A xf / 

Explicit Representation (n = 2 to N~ 1) (4) 

where the superscript i refers to the current iteration at time t and the 
superscript /+ 1 refers to the next iteration at time t + At. It can be seen that 
with the explicit representation, C r „ + 1 can be calculated directly, or explicitly, 
from the known concentrations at time /, C' n . A portion of a concentration 
grid at time t and the new concentrations at t + At, are shown schematically 
in Fig. 1 . Forward and backward difference equations, such as for calculating 
the flux at a boundary (i.e., at n= 1 or n = N), and higher order difference 
equations containing more terms than those shown in Eqs. (3) and (4) (e.g., 
C„- 2 and C fi+2 h are available. 910 However, higher-order difference equations 
have seldom been used in modeling diffusional transport associated with 
high-temperature corrosion. The time increment for the explicit method is 
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Fig. 1 . Schematic concentration grid for time i (squares) and time t + At 
(circles). The concentration subscript refers to the node number and the 
superscript refers to the time iteration. 


limited by a stability criterion 9 given as 
A t 

D -—^<0.25 Stability Criterion for Explicit Method (5) 

The most common implicit F-D method incorporates the Crank- 
Nicholson (C-N) representation. For the case where Z)//(C), Eq. (2b) 
becomes 

= i D } (c;V,-2c^ l + c; 1 t 1 1 ) { (c;, +1 - 2 c;,+c;, 1 ) { 

A t 2 I A.Y 2 AX 2 j 

Implicit C-N Representation (n = 2 to N) (6) 

It can be seen that with the implicit representation, calculating CJ/ 1 requires 
not only the known concentration values at time /, but also the unknown 
values C‘ t ^\ and CJ,Vi . However, it can be shown 9 that when Eq. (6) is used, 
a very desirable tridiagonal matrix of unknowns can be formed which can 
be solved using standard matrix algorithms to determine the values for the 
CJ, + 1 . The C-N implicit method is very stable and does not suffer from the 
restriction on the size of At as does the explicit method. Larger time incre- 
ments are desirable since they result in less iterations, less computation time 
and less possibility of roundoff or truncation error. 

When diffusion coefficients are strongly concentration dependent, each 
of the differential terms in Eq. (2c) are replaced by F-D equivalents and an 
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explicit method of solution must be used subject to the stability criterion 
for At given in Eq. (5). An explicit method may also be desirable when a 
strongly time-dependent boundary condition requires the use of small time 
increments negating some of the benefit of the greater At of the implicit 
method. In addition, the explicit method is also attractive by generally being 
more straightforward to program and routinely yields satisfactory results 
when computer time, or the likelihood of significant roundoff error, is not 
a concern. Otherwise, whenever the concentration dependence of the diffu- 
sivities is unknown or can be ignored, the implicit F-D method is to be 
preferred since it is inherently more stable and thus a more efficient solution 
allowing the use of larger values of At. 


Initial and Boundary Conditions 

A specific solution to Fick’s second law, whether by explicit or implicit 
methods, requires initial and boundary conditions. Concentration values are 
initially assigned to each of the nodes in the grid to reflect the starting 
concentration condition within the sample or component. This initial condi- 
tion may be as simple as a uniform concentration throughout a sample (e.g., 
C„ — C for n- 1 to N. w here C is the bulk concentration of the component ), 
or a complicated profile, such as a near-surface solute enrichment due to 
the presence of a coating 11 ' 12 or a near-surface solute depletion resulting 
from an intentional, or unintentional, preoxidation exposure. 13 14 

Boundary conditions are required to calculate the new boundary con- 
centrations, C \ + 1 and C!v \ These boundary conditions are generally defined 
by either a fixed concentration or a flux condition. A constant solute concen- 
tration at the oxide metal interface during isothermal oxidation (e.g., C\ = 
C s ) is an example of a fixed concentration boundary condition at an outer 
surface. A constant solute concentration at the innermost node which never 
changes due to the thickness of the sample or component (e.g., C N =C ) is 
an example of the same type of boundary condition applied at an inner 
boundary. This latter boundary condition is appropriate for samples which 
can be considered "infinitely thick’' in comparison to the diffusion distances 
during the time span of the corrosive attack. 1 ^ In contrast, diffusion often 
reaches to the center of a sample or component. In this case (i.e., finite 
sample thickness), a zero-flux plane is commonly used as an inner flux 
boundary condition. 16 A zero-flux plane indicates that solute is being 
removed from, or corrodant is being supplied to the center of the sample in 
both directions equally. Similarly, a flux of N or C into an alloy controlled 
by the dissociation of ammonia or methane on the surface is an example of 
a flux boundary condition at an outer surface. 1 ' ls 
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Moving Boundaries, Fixed- and Flexible-Grid Schemes 

Many problems associated with high-temperature corrosion processes 
involve moving boundaries. Examples of such moving boundaries include 
recession or growth of an oxide metal interface, growth of a subscale layer, 
and growth of a different, near-surface phase in the alloy. These moving 
boundaries have been handled in different ways by F-D studies. In many 
problems, a flexible-grid scheme is used in which the total grid expands or 
contracts to match the interface motion. When the grid expands or contracts, 
the concentrations at the nodes must be adjusted. The most common method 
of shifting the node positions and adjusting the node concentrations is that 
referred to as a “Murray-Landis” (M-L) transformation .' 9 This transforma- 
tion shifts the nodes an amount proportional to their position from the 
moving boundary to maintain a uniform node spacing. Although the M-L 
transformation has been used extensively in the past, Crusius et al 2i) suggest 
an extensive modification which yields a higher accuracy under a wider range 
of conditions. However, this modification results in a considerable increase 
in the computation time. Moving boundaries using a flexible-grid scheme 
and an M-L transformation can be treated by either implicit or explicit F- 
D methods. Several examples in the applications section of this paper have 
employed the M-L transformation . 1 1,1215 16 In contrast, one study has also 
shown the feasibility of modeling interface motion with a fixed-grid scheme . 21 

Lastly, the formation and growth of carbide subscale layers has been 
handled in a different fashion by considering the carbide formation at each 
individual node. Two methods have been used, both with fixed grid schemes. 
These methods are discussed in greater detail in the application section 
dealing with carbide formation. 

APPLICATIONS OF THE F-D TECHNIQUE 

The following sections are intended to briefly review several applications 
of the F-D technique in the area of high-temperature corrosion with the 
intent to give the reader a sense of the applicability and flexibility of the 
technique. Two general cases will be examined, firstly, corrosion processes 
that involve solute transport outward, such as occurs in the selective oxida- 
tion of Al- or Cr-containing alloys and coatings, and secondly, corrodant 
diffusion inward, such as occurs during nitridation, internal oxidation and 
carburization. 


Case I. Solute Transport 

Protective External Scale Formation 

Whittle et al . 22 presented one of the earliest known studies utilizing the 
F-D technique to predict solute transport in binary alloys. Their model 
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examined the effect of parabolic, cubic and logarithmic oxide growth rates 
on the solute concentration and time dependence of the solute concentration 
at the alloy surface, C b Hence, the flux boundary condition in terms of the 
solute B can be stated as 




= —D ^ 

X = 0 cX 


- rj(kit + k u ) ^ 


(7) 


where J B is the solute flux in the alloy at the oxide-metal interface, k j and 
k u are oxidation rate constants, rj is a constant dependent on the oxide 
stoichiometry, the units for the solute flux and oxide growth rates, etc., and 
<p has the value \ for a parabolic rate, f for a cubic rate, and 1 for a 
logarithmic rate. The model quantified the time dependence of the solute 
concentration at the alloy surface and predicted concentration-distance pro- 
files within the alloy for the three oxide growth rates. In applying the model 
to the isothermal oxidation of various Fe-Cr alloys, reasonable agreement 
was observed between measured and predicted Cr concentrations at the 
oxide metal interface and concentration-distance profiles within the alloys 
at different times. 23 This analysis was then applied to show that the break- 
down in protective Cr^O^ scale formation on Fe-(13- 25 wt.%)Cr alloys can- 
not be attributed to a chemical transformation within the alloy due to a 
general depletion of Cr since the Cr concentration is never reduced below 
that necessary for thermodynamic stability of the Cr 2 0^ scale. 

Nesbitt et al. also used F-D techniques to predict A1 concentrations in 
P NiAl alloys, 16 and Al and Cr concentration profiles in y NiCrAl alloys 15 
undergoing cyclic oxidation. During thermal cycling, the oxide which grows 
at high temperature undergoes partial spalling on cooling. The practical 
consequence of this oxide spalling is to increase the rate of solute consump- 
tion over that for the isothermal case. An oxide spalling model, COSP, 24 
which predicted this accelerated loss of Al, was used as a flux boundary 
condition at the oxide metal interface for both the NiAl and NiCrAl studies. 
In addition, both studies used an explicit F-D method to take into account 
the concentration dependence of the diffusivities (see Eq. (2c)). Recession 
of the alloy surface due to the loss of Al during selective oxidation was also 
taken into account. This oxide- metal interface recession, A£, is defined as 15 

A£ = -K ai 72|A/ (8) 

where V A \ is the partial molar volume of Al in the alloy and Jm is the flux 
of Al entering the oxide (which is greater than the flux of Al in the alloy at 
the oxide metal interface because of the interface motion). An M-L trans- 
formation was used to adjust the concentration grid for this alloy surface 


recession. 



Numerical Modeling of High-Temperature Corrosion Processes 


317 


In the NiAl modeling study, 16 the F-D model predicted A1 concentration 
profiles within thin coupon samples after various cyclic oxidation exposures. 
A zero-flux boundary condition was used at the center of the sample where 
the A1 concentration was shown to decrease due to the limited finite sample 
thickness. The modeling showed that the A1 concentration throughout the 
entire sample (up to 2 mm thickness) decreased somewhat uniformly as A1 
was removed from the sample by oxidation. These predictions were verified 
by concentration measurements on oxidized samples. This relatively uniform 
decrease in A1 was due to the rapid A1 transport in the alloy in comparison 
to the rate that A1 was consumed by oxidation. The practical consequences 
of this observation are (1) oxidation cannot be viewed as only affecting a 
near-surface region, and (2), the mechanical properties of P NiAl alloys 
which are concentration dependent can be significantly affected by the 
uniform loss of A1 during oxidation. Cyclic oxidation testing of NiAl alloys 16 
showed that the protective ALCL scale was slowly replaced by a less- 
protective NiAl 2 0 4 scale when the A1 concentration in the alloy decreased 
below approximately 36 at.%. This value was then taken as a failure criterion 
and the protective lifetimes of various p NiAl alloys of different thickness 
were predicted. Reasonable agreement was shown between predicted and 
observed protective lifetimes. 

In NiCr A1 alloys, 1 5 both A1 and Cr concentration profiles were predicted 
for the case of cyclic oxidation. For NiCrAl alloys, Fick’s laws can be written 
as 


JJ vX J ' dX 


j , k = Al, Cr 


First Law (9) 


and 


PCj = d( Djj(dCj/vX)) + d(Dj, k ( dCk/dX )) 
dt dX dX 


j\ k - Al, Cr 


Second Law (10) 


where the main diffusivity, D ) h relates the concentration gradient of compo- 
nent / to the flux of j and the off-diagonal, or cross-term diffusivity, D jk , 
relates the concentration gradient of component k to the flux of j. Equation 
(10) was expanded as for the binary case (Eq. (2c)) and each of the terms 
was replaced with a F-D equivalent for solution by the explicit method. The 
oxide-metal interface recession due to Al loss was calculated as shown above 
for the NiAl alloys (Eq. (8)). However, since Cr is not consumed in the 
oxide scale, Cr must diffuse away from the receding oxide-metal interface 
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(X = c , ). The outer boundary condition for Cr then becomes 

( 11 ) 

(It 

where J Cr is the flux of Cr in the alloy at the interface and Co is the Cr 
concentration at this interface. 

The F-D model 15 predicted the rate at which the A1 concentration at 
the oxide metal interface decreased with time due to the higher rate of A1 
consumption associated with cyclic oxidation. The model was also used to 
predict the minimum A1 necessary to maintain a protective Al 2 Ch scale 
during cyclic oxidation. It is well known that the addition of Cr to NiAl 
alloys greatly reduces the amount of A1 necessary for protective scale forma- 
tion during cyclic oxidation. Good agreement was observed between pre- 
dicted and experimental data for alloys with Cr concentrations greater than 
20 at.%Cr. However, the model significantly underpredicted the amount of 
A1 for binary Ni A1 alloys by a factor of 2.5. The modeling results showed 
that the beneficial effect of Cr additions to Ni A1 alloys could not be attri- 
buted to any transport-related effect of Cr enhancing the A1 transport to 
the oxide scale. The implication of this result is that the addition of Cr likely 
acts as a sacrificial element reducing the transient formation of Ni oxides 
when the bare metal is exposed, either at the start of oxidation, or after 
oxide spallation. 


Overlay Coating Degradation 

Degradation of y + fi MCrAl (M = Ni or Co) overlay coatings on Ni- 
and Co-based substrates during isothermal and cyclic oxidation has also 
been modeled in two studies. 11 12 Overlay coatings are deposited on super- 
alloys to provide an A1 reservoir and extend the protection afforded by the 
protective Al 2 Cb scale. During high-temperature exposure in an oxidizing 
environment, this reservoir of A1 in the coating is reduced both by consump- 
tion of A1 by selective oxidation and by interdiffusion of A1 into the sub- 
strate. Therefore, A1 diffusion to the oxide and into the substrate was 
modeled. The two models used different failure criteria to quantify the pro- 
tective life of the coating. This selection of the failure criterion affected how 
diffusion within the y + /? coating was modeled. 

In the earlier model by Nesbitt and Heckel 11 the protective coating 
lifetime was defined as the time at which the A1 concentration at the oxide 
coating interface decreased to zero. Experimental studies showed that the fi 
phase within the coating was totally depleted long before the A1 concentra- 
tion at the oxide metal interface approached zero. Consequently, it was 
assumed that the two-phase structure of the coating could be represented 
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Fig. 2. Schematic Al concentration profile after oxidation of a 
Y + P MCrAlY overlay coaling on a superalloy substrate. The 
dashed line represents predictions from the model by Nesbitt cl 
al and the solid line represents predictions for the model by Lee 
cl al. ]Z Arrows indicate the direction of Al diffusion from the y + /? 
region of the coating. 


as a single /-phase layer with the initial condition given as the initial Al and 
Cr concentrations of the coating (Cm, C tr ). The Al concentration profile 
in the coating and substrate showing the effect of this assumption is shown 
schematically in Fig. 2. Since diffusion of Al to the oxide scale and into the 
substrate is largely controlled by diffusion in the / phase, concentration- 
dependent, ternary dilfusivities appropriate for the / phase were used. Cyclic 
oxidation conditions were modeled by taking the rate of Al consumption 
predicted by the COSP oxide spalling model 24 as a flux boundary condition 
at the oxide coating interface. Surface recession due to loss of Al to the 
oxide scale was also taken into account. Good agreement between measured 
and predicted concentration profiles was observed after both cyclic furnace 
testing of NiCrAlZr coatings on simple NiCrAl substrates and cyclic burner 
rig testing of NiCoCrAlY coatings on a Ni-base superalloy (Fig. 3). 2:1 Based 
on this satisfactory verification of the model's predictive capability, the effect 
of various system parameters (e.g., the initial Al and Cr content in the 
coating and substrate) on the protective life of the coating was examined. 
A further practical consequence of this study was the demonstration that 
more Al can diffuse into low Al-containing substrates than is lost to the 
oxidation process. 

In the later model developed by Lee el al.J 2 the protective coating 
lifetime was defined as the time at which the /? phase was totally depleted 
from the coating. The practical application of this failure criterion might 
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Fig. 3* Measured and predicted A1 concentration profiles 
for a NiCoCrAlY coating on U-700 superalloy after 1000 
1-h cycles in a furnace at 1100 CY 


signal the time at which a coated component be removed from service and 
the coating stripped and a new coating applied. Recession of the /? phase 
due to A1 diffusion from the coating to the oxide scale and into the substrate 
was accomplished by establishing moving boundaries on both sides of the 
coating (Fig. 2). A flexible-grid scheme was used and an M-L transformation 
was employed to shift the node positions and concentrations after the inter- 
face motion. The boundary condition at the y/y + P interface is given by 

( Cai - Ctf + p ) d$/dt = JW ' I x - ; (1 2) 

where J r \\ 7 * p is the flux in the y phase at the y/y- f /? interface located at 
X = £, Cai is the initial A1 concentration in the coating, and C is the 
maximum A1 concentration in the y phase given by the endpoint of a tie 
line passing through CXi in the y + fi region of the phase diagram. No 
diffusion was allowed in the y + fi region of the coating and the value of 
Cai 7 ^ on each side of the coating was assumed constant. The A1 concentra- 
tion dependence of Aaiai was taken into account although the effect of Cr 
on the diffusion of A1 in the y layer was ignored (i.e., only A1 transport was 
modeled). A rate of A1 consumption appropriate for either isothermal or 
cyclic oxidation was used as a flux boundary condition at the oxide- coating 
interface. 24 Good agreement was observed between the measured and pre- 
dicted volume fraction of P remaining in the y + P region of the coating 
following isothermal oxidation of a CoCrAlY coating on a Co-base super- 
alloy (Fig. 4). This model was also used to evaluate the effect of various 
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Fig. 4. Measured and predicted volume fraction of 
P phase in the y + p Co 17Cr 23A1 0.3Y (at.%) 
coating on two superalloy substrates after isother- 
mal oxidation. 12 


parameters (coating thickness, oxide growth rate, isothermal vs. cyclic oxida- 
tion conditions, etc.) on the protective life of the coating for the failure 
criterion stated above. It is noteworthy that both models simulating overlay 
coating degradation were used to perform parametric analyses to examine 
the impact of various parameters on the coating life. This use of these models 
emphasizes an important aspect of modeling studies, namely, to use the 
model to identify the most important parameters having the greatest impact 
in a problem, especially in complex problems such as multiphase coatings 
on superalloy substrates. 

Case II. A Corrodant Transport in the Absence of Internal Precipitation 

Carburization, or nitridation, of Fe- and Ni-based alloys may be 
intentional, such as to provide some near-surface mechanical properties in 
a component, or unintentional, such as occurs upon exposure of components 
during high-temperature operation in a C- or N-containing environment 
(i.e., a high-temperature corrosion process). Modeling the carburization or 
nitridation process, whether it occurs intentionally or as the result of corro- 
sion, is accomplished in an identical manner. Hence, several examples below 
involve the intentional nitridation and carburization of Fe-based alloys and 
are discussed to demonstrate the flexibility of the F-D technique. 

Nitridation 

Nitriding a (bcc) Fe-base alloys is a common practice to improve the 
tribological properties of the surface by the formation of a near-surface. 
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nitride-containing layer and an inner region containing N dissolved interstiti- 
ally. Analytical solutions to this diffusion problem exist if the surface concen- 
tration, Cn, is assumed to be constant. However, Rozendaal et ai xl 
presented experimental evidence that Cn only slowly approaches the equilib- 
rium concentration during nitridation. Hence, these authors used F-D tech- 
niques to model the diffusion of N into thin sheets of u-Fe. The time to reach 
a critical surface concentration where y f Fe 4 N would form was predicted, but 
the precipitation and growth of this nitride subscale layer was not modeled. 
The nitriding process based on the dissociation of ammonia, can be broken 
into two steps (assuming no lack of ammonia transport to the surface): (1) 
the dissociation of ammonia at the Fe surface, and (2) diffusion of N into 
the sample. The dissociation process is reaction-rate controlled and can be 
described as 


Js — A ' ( C n 1 ( n) 


(13) 


where Js is the flux of N, A' is the reaction rate coefficient, C s is the N 
concentration in ar-Fe in equilibrium with the gas atmosphere, and Cn is 
the actual N surface concentration. Hence, the boundary condition at the 
surface can be derived by equating this flux of available N on the surface 
with the flux of N diffusing into the alloy given by Fick's first law, such that 


-D 


rC s 

cX 


= k{Cs - C\ ) 


(14) 


Rozendaal et ai xl employed the implicit C-N representation (Eq. (6)) to 
predict N concentration profiles after various nitriding exposures. Good 
qualitative agreement was shown between predicted concentration profiles 
and microhardness profiles indicating a time dependence of the surface con- 
centration Cn. The predicted time dependence for various ammonia mix- 
tures, and specifically the time to reach a critical surface concentration when 
y Fe 4 N would form, was used to explain the incubation time experimentally 
observed prior to y' Fe 4 N formation. 


Oxidation ami Oxygen Dissolution 

Denis et a/. 21 developed an implicit F-D model to simulate O transport 
through an external scale and into the matrix. This model, revised from an 
earlier explicit model by Garcia 26 in order to decrease computation times, 
presents a rare example of the use of a fixed-grid scheme with a moving 
boundary. Concentration-independent diffusion coefficients were assumed. 
The model was applied to the oxidation of Zircaloy-4 in water vapor which, 
as with many Zr alloys, initially obeys a parabolic law. However, continued 
oxidation results in a sharp transition to linear kinetics which are attributed 
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to physical changes in the oxide (Zr0 2 v F such as pore and crack formation. 
The time interval each iteration. At, was chosen by an iterative process such 
that the new r oxide metal interface position each iteration corresponded to 
a node location (i.e., a time increment was chosen such that the interface 
motion. As, was always equal to the node spacing, AX ). In brief, this itera- 
tion process consisted of calculating new concentrations, Cc^ 1 , using Eq. 
(6), and interlace velocities, d£/dt, using a flux balance equation similar to 
Eq. (12), for different values of At. An equation was derived relating the 
interface velocities and time increments w hich allowed an estimate of a new r 
At such that Ac^AA\ This process w'as repeated until a value of At was 
achieved which produced an interface motion essentially equal to the node 
spacing. 

For convenience, O transport during the parabolic portion of the oxida- 
tion of Zircaloy w r as performed using analytical equations. Modeling the 
linear kinetics w r as accomplished in two w'ays. First, a time-dependent value 
for Do, the diffusion coefficient in the zirconia scale, was used. To determine 
this time-dependence of D 0 , the model was operated in a manner that a 
value for D () was determined each iteration such that the predicted sample 
weight gain, which was proportional to the predicted oxide thickness, 
matched that previously measured for the oxidation of Zircaloy. At a given 
temperature, the derived values for D () , which were constant in the parabolic 
region, instantaneously increased at the transition from parabolic to linear 
kinetics, and increased linearly with time thereafter. Although time- 
dependent diffusion coefficients have little physical significance, this example 
illustrates the flexibility of the F-D technique to handle time-dependent 
parameters. 

The linear kinetics were also modeled by assuming a two-part oxide 
scale. 21 The model was revised to consider diffusional transport through an 
inner region of the Zr0 2 scale only, while the outer region was assumed to 
offer no resistance to O transport due to formation and growth of pores 
and channels. The diffusion coefficient in the oxide was constant during both 
the parabolic and linear oxidation regimes. At the time of the parabolic/ 
linear transition, an interface was introduced in the oxide as the boundary 
between the porous layer containing channels and the dense inner layer 
which controlled diffusion, A schematic concentration profile indicating 
these layers is shown in Fig. 5. It was also assumed that the O potential 
within the channels was the same as that of the environment such that the 
O concentration throughout the porous oxide layer was assumed constant 
and equal to the concentration at the gas oxide interface, Co. Since the 
weight of the sample increased linearly, and therefore the flux of O into the 
oxide was constant, the width of the dense inner oxide layer remained con- 
stant and moved into the Zircaloy at a constant rate. The usefulness of the 
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Fig. 5. Schematic O concentration profile showing the two- 
layer Zr0 2 v scale on Zircaloy used by Denis cl al: ] 


model was in determining the value for Z) 0 during the parabolic oxidation 
regime, and in developing a plausible model based on a dense inner oxide 
layer of constant thickness to explain the linear oxidation regime. 


Carburization 

Goldstein and Moren 13 presented one of the earliest studies employing 
F-D techniques to model C diffusion in Fe C and Fe C X alloys, where X 
is a solute element such as Ni, Mn, Cr or Si. Each of these elements is 
known to affect both the diffusivity and solubility of C in austenite. For the 
binary case, various carburization treatments were modeled where both the 
temperature and surface C concentration, Q , were time dependent. The 
temperature and concentration dependence of the diffusion coefficient of C 
in austenite, /)£, were used in the model. Flowever, since Dl is not a strong 
function of concentration, the term dD/vC(dC/dX ) in Eq. (2c) was ignored 
so that concentration dependent values for D( were used with Fick’s second 
law given by Eq. (2b). A fixed-grid scheme was used with a constant number 
of nodes. Predicted and measured concentration profiles are shown in Fig. 
6a. The variation of Cl and temperature with time are shown in the inset. 
The agreement between the measured and predicted C concentration is 
apparent. 

For the ternary alloys, Goldstein and Moren 13 gave the solute concen- 
tration (Cr or Si) a steep concentration gradient near the surface, such as 
might occur due to limited oxidation prior to carburization, to evaluate the 
effect of the solute gradient on the diffusion of C. In their work, the C 
content of the steel was fixed at 0.15 wt.% with a surface concentration of 
1.0wt.%. The solute concentration, either Cr or Si, was set at 2 wt.% or 
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Fig. 6. (a) Predicted and measured C concentration profiles for the 
carburization treatment shown in the inset, and (b) predicted C con- 
centration profiles for Fe 2Cr 0.3 5C, Fe 2.5Si O.lSCandFe 0.1 5C 
(wt.%) alloys carburized at 1700 F for 8.31 h. Surface C concentra- 
tion Cl = 1.0. 13 


2.5 wt.%, respectively, throughout the sample but zero at the surface node. 
Chromium, a strong carbide former, has a negative cross-term diffusion 
coefficient which is approximately one third the value of the main diffusion 
coefficient, or A:,c r /A:\t = — 0.34 (i.e., a Cr concentration gradient opposite 
in sign to that of C will increase the flux of C [see Eq. (9)]). This relatively 





326 


Nesbitt 


strong cross-term effect, coupled with a large, positive Cr concentration 
gradient near the surface, resulted in a maximum in the C concentration of 
nearly 1 .7 wt.%. The practical effect of this maximum is that more C diffusion 
occurs into this alloy than that for the binary alloy. In contrast, the cross- 
term effect in Si-containing alloys is positive (i.e., 0 ( Sl /D r . r = +0.30) such 
that the large positive Si concentration gradient near the surface caused a 
sharp decrease in the near-surface C concentration (i.e., the Si concentration 
gradient opposes the influx of C). The predicted C concentration profiles 
for these three alloys, Fe C, Fe C Cr and Fe C Si are shown in Fig. 6b. 

Case II. B Corrodant Transport Accompanied by Internal Precipitation 

The inward transport of a corrodant often results in internal precipita- 
tion and subscale formation due to reaction of one or more alloying elements 
with the corrodant. The manner in which this precipitation is handled by 
F-D techniques depends, in large part, on the thermodynamic stability of the 
reaction product. Internal oxidation presents an example of the formation of 
a compound with high thermodynamic stability. Carbide formation during 
carburization presents an example of the formation of a compound with 
moderate thermodynamic stability. 

Several of the following studies involving internal precipitation make 
use of a “labyrinth factor” to account for the reduction in cross-sectional 
area due to precipitate formation. This factor, which has a value between 
zero and one, is used to reduce the matrix diffusivity. The value of the 
labyrinth factor is generally related to the volume fraction of the remaining 
matrix. 


Internal Oxidation 

Internal oxidation involves the diffusion of O into an alloy and the 
formation of a subscale of solute oxide precipitates (BO v ). The high stability 
of many oxides results in an extremely low solubility product such that the 
concentration of solute in solution within the subscale and the O concentra- 
tion at the subscale alloy interface (i.e., at the reaction front) is nearly zero. 
Furthermore, little O penetrates beyond the precipitation zone into the alloy. 
Hence, for a highly stable oxide, it is generally assumed that (1) the oxide 
precipitation reaction occurs only at the subscale alloy interface (i.e., at X — 
£); (2) all the solute in the subscale is converted to oxide (i.e., C B = 0 within 
the subscale); and (3) the concentration of solute and O at the subscale 
alloy interface is zero (i.e., C B = C o = 0 at X= £ ). A grid is established for 
O in the subscale and, if B is sufficiently mobile, for B in the alloy. (The 
diffusivity of alloying elements can be several orders of magnitude less than 
that for corrodants which diffuse interstitially (e.g., O, C and N) such that 
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Fig. 7. Schematic concentration profiles for O (G>). 
solute B (Oh), and B precipitated as BO, (C{l p1 ) for the 
internal oxidation of a B-containing alloy where B forms 
a highly stable oxide BO, . 7 


the alloying elements are often assumed to be immobile). Schematic concen- 
tration profiles for the outward diffusion of B and the inward diffusion of 
O are shown in Fig. 7. The flux of O and the flux of B arriving at the 
subscale alloy interface are related according to the stoichiometry of the 
precipitated oxide BO, such that 


D B 


SC H 

dX 


Do c C 0 
v (X 


(15) 


where D 0 is the diffusion coefficient of O in the alloy and v is the stoichiome- 
try factor. The interface motion in the time interval At is given by a 
simple mass balance relationship which can be stated as 


(~ppt 4.1 = _ P.v 

At v d X q v 


(16) 


where Cfi pt is the concentration of precipitated solute in the subscale (Fig. 
7). The O profile through the subscale is often assumed to be linear, which 
allows the O concentration derivative term dC a to be replaced with the O 
difference across the subscale (i.e., Co 0) and the distance derivative oX, to 
be replaced w'ith the w'idth of the subscale, £, yielding Co/£ for the O 
concentration gradient in Eq. 16. If B is immobile, the concentration of B 
in the subscale CE pt (existing as BO,.) remains equal to the concentration of 
B in the alloy Cr. However, if reasonable diffusion of B does occur, Cr p ‘ 
in the subscale will exceed Cr in the alloy by an enrichment factor typically 
defined as a = Cr^/Cr. 
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Vedula et al 21 developed a F-D model to describe internal oxidation 
for planar, cylindrical and spherical geometries of finite size where the solute 
B was considered mobile. The model was initially used to perform parametric 
analyses examining the subscale thickness and enrichment under isothermal 
conditions with a constant boundary condition (Cq constant), and was later 
extended to consider the effect of nonisothermal conditions and a time- 
varying boundary condition (Co =/(/)). The concentration profiles in the 
alloy and subscale were determined by an explicit F-D method. Each itera- 
tion consisted of ( 1 ) calculating new concentrations at each node in the 
subscale and alloy according to Fick’s second law, (2) calculating the inter- 
face motion A£, (3) shifting the nodes in both the subscale and alloy and 
correcting the concentrations according to an M-L transformation, (4) calcu- 
lating Cb p1 in the subscale (for a > 1 ). 

The interface velocity d^/dt is directly related to the O flux into the 
alloy which is directly related to Co. It had previously been observed by 
Rapp 3 that the interface £ reversed directions when Cq was quickly 
decreased. A second, more gradual interface reversal eventually occurs caus- 
ing the interface to again move in the original direction into the alloy. This 
double interface reversal produces a localized band of high enrichment, 
known as an interruption band. This interface reversal was also quantita- 
tively modeled by Vedula et al 21 for a sheet of thickness L. The normalized 
interface position, displaying an interface reversal of A £/L, is shown in Fig. 8 
(bottom). The rate of precipitate formation is directly related to the interface 
velocity. After the second reversal, the interface again moves into the alloy 
at a slower rate in accordance with the lower O flux. The reduced interface 
velocity (d^/dt ) 2 results in a higher rate of oxide precipitation, as shown in 
Fig. 8 (top). The second, more gradual interface reversal results in a “spike” 
in the rate of oxide precipitation. 

Carburization Accompanied by Carbide Formation 

The moderate stability of many carbides results in a much greater solu- 
bility product for the carbides in comparison to the oxides. Consequently, 
there can be significant C diffusion into the alloy beyond the carbide subscale 
and a significant concentration of the carbide-forming element within the 
subscale region. Since in many alloys of interest several carbides of the same 
element may form, it is useful to consider a carbide stability map such as 
that shown schematically in Fig. 9a for a simple binary alloy which forms 
carbides I, II and III (e.g., a Ni-Cr alloy with carbides Cr 23Q, Cr 7 C} and 
O3C2). Initially, consider a low C activity in the alloy, such as showm as 
point P in the figure. During carburization, C will diffuse into the alloy 
raising the activity of C, as indicated by the dashed line PQ. At point Q, 
carbide I begins to precipitate internally starting at the surface. Further 
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Fig. 8. Time versus distance plot showing a double inter- 
face reversal (bottom) and the associated rate of oxide 
precipitation (top) during internal oxidation. The inter- 
face velocities, where (dg/dt) ] > (rfc /dt) 2 , are also 
shown. The arrows indicate the path direction . 27 

diffusion of C into the alloy will result in the inward growth of this carbide 
subscale layer while carbide I continues to form and grow within the layer 
itself (unlike that for internal oxidation where no further oxidation occurs 
in the subscale layer). If the C concentration at the outer surface remains 
constant, the subscale layer will thicken parabolically with time. As further 
carbide l formation occurs in the subscale, the activity and concentration 
of B in solution in the subscale decreases along the line indicated as QR in 
Fig. 9a. When the C activity at the outer surface reaches that indicated by 
point R, a carbide II layer will begin to form and grow inward as carbide I 
is converted to carbide II at the carbide I /II interface. Continued diffusion 
of C into the alloy will cause the continued growth of both carbide subscale 
layers along the lines QR (inner layer of carbide I) and along RS (outer 
layer of carbide II). If the C activity outside the alloy is sufficiently high, 
such as for dc shown in Fig. 9a, a third carbide (III) may form at the outer 
surface and move inwardly with the C activity in subscale III following the 
line ST. 
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Fig. 9. C'a rburization of a binary alloy where three carbides 
may form with a single element B (a) schematic carbide 
stability map and (b) schematic concentration profiles where 
subscales of three carbides in (a) are shown. ( After Sockel el 
al. 2 *) 

Schematic concentration profiles for the carburization of a binary 
alloy, based on the modeling predictions for Ni 25Cr by Sockel et al , 28 are 
shown in Fig. 9b. The solute B was assumed immobile. There is a significant 
concentration of the solute B within the subscale and a small but apparent 
amount of C in the alloy below the carbide I subscale resulting from the 
moderate stability of the carbides. The C activity, indicating the activities 
of interest from Fig. 9a (P, Q, R, S, T), is also shown schematically in 
Fig. 9b. 

Two methods appear in the literature to model, by F-D techniques, 
the carburization process when accompanied by carbide formation. Both 
methods handle formation of multiple carbides of the solute B. The first 
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and earlier method, used by Bongartz et r//., IH uses a two-step procedure at 
each node to calculate the concentration of B and C in solution and in the 
carbide for the current time iteration At. The two steps involve, firstly, a 
diffusion calculation, and secondly, an equilibrium calculation to partition 
the solute and C between that which remains in solution and that which 
forms the carbide. The first step calculates the C diffusion to a node in the 
time interval At according to Fick’s second law. Typically, this step increases 
the C concentration of the node. The second step compares the product of 
the current C and B concentrations with the solubility product for the car- 
bide, K s . The excess C, AC, above that given by K s , is assumed to form 
carbide. The concentration of soluble C at the node is reduced by AC while 
the concentration of C in the carbide is increased by AC. The quantity of 
B w hich is involved in the carbide formation with AC is also removed from 
the concentration of soluble B at the node. The new C concentration of the 
node, which now satisfies K s , is assumed available to diffuse to the next 
node where this two-step process is repeated. The activities associated w'ith 
these concentrations of C and B at this particular node will correspond to 
a point on the line QR on the carbide stability map shown in Fig. 9a. This 
two-step procedure begins at the outer nodes and progresses inward. As 
time progresses and the concentration of B in solution decreases w r hile the 
concentration of C increases, the activities will follow the line from Q to R. 
When point R is reached, carbide II will also be stable at this node. Further 
diffusion of C to this node will only result in the transformation of carbide 
I to carbide II without a further change in the concentrations of B or C 
until all of carbide I at the node has been transformed. Thereafter, further 
diffusion of C to this node will cause additional formation of carbide II 
while the concentrations, and activities of B and C in solution at the node 
now' follow along the line RS in Fig. 9a 

This two-step procedure was extended to the problem where many 
carbides of different solute elements may exist simultaneously, as for the 
case of the carburization of various Fe Ni Co Cr Mo alloys. 29 For the 
second step in this case, the excess C, AC, which exceeds the solubility limit 
of the various possible carbides must be partitioned between those carbides. 
For M possible carbides, the excess C can be described as AC = £AC M w'here 
AC m reflects the C tied up in each of the various carbides which can contain 
varying amounts of M, w'here M stands for the carbide-forming element 
(e.g., Cr, W, Ti). By assuming a regular solution model for the matrix, a 
series of M nonlinear equations, based on concentrations known at time t u 
can be derived to determine the M values of AC M . Numerical methods were 
used to solve for each of the AC M ’s. This procedure is used repetitively at 
each node to determine the C in solution and that in each of the carbides 
as a function of position in the sample. 
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Fig. 10. Measured and predicted C concentration 
profiles following carburization of a Ni lOCr 
alloy. The various predicted curves represent the 
indicated values for the solubility product K , for 
CnC, and D C M) 


Bongartz et al. M) applied these models to the carburization of various 
Ni-Cr alloys. Although the necessary thermodynamic and kinetic data 
(primarily the solubility products for each carbide, activity coefficients, y c , 
y Cr , and c ) exist in the literature, these published data exhibit significant 
scatter, ft was shown that the calculated C profiles are very sensitive to the 
thermodynamic and kinetic data, as shown for a Ni lOCr alloy in Fig. 10. 
Hence, Bongartz et a /. matched predicted and measured C profiles to get 
the most reasonable value for D cx and the solubility products for the Ni 
lOCr alloy. The values for D C \ C , 7c and the solubility product for Cr 3 C 2 
determined from the Ni lOCr alloys were then applied to gain the solubility 
products for Cr 7 C 3 and Cr 2 3 C 6 . Because of the sensitivity of the calculated 
profiles to the input data, the authors clearly demonstrated the usefulness of 
the F-D technique, when combined with reliable experimental concentration 
measurements, to derive important thermodynamic and kinetic data. 

Bongartz et al . ] 8 also simulated carburization and carbide precipitation 
in Fe-based alloys with planar, cylindrical or spherical geometries. A con- 
stant C flux, determined as an average value from the C pickup on experi- 
mental samples, was taken as the outer surface boundary condition. Similar 
to the systematic approach stated above, predicted and measured C concen- 
tration profiles in Fe Ni Cr alloys were used to derive an approximate 
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diffusion coefficient for C in this alloy. Likewise, the solubility product for 
NbC was then derived by comparing measured and predicted total C concen- 
tration profiles (C in solution and as carbides) in Fe-Ni-Cr-Nb alloys. The 
model was then used to make semiquantitative predictions of the degrada- 
tion of Fe Ni Cr Nb tubes exposed to a carburizing environment at 900 C 
for times up to 100,000 h (i.e., 15-year operation in a petrochemical plant 
or in the steam reformer of a nuclear process heat system). The predictions 
were able to quantify the difference in carburization attack for the condition 
of a constant C concentration at the outer surface (i.e., a high and unlimited 
supply of C from the environment controlled only by diffusion into the 
alloy) versus a constant flux of C at the outer surface (i.e., the more common 
situation where the rate of supply of C to the alloy is limited by the decompo- 
sition of the carburizing species (e.g., methane) at the alloy surface or gas- 
oxide interface if an oxide scale is present). 

Bongartz et al , 29 later utilized the same systematic approach to model 
the carburization of various binary, ternary and quaternary Ni-Fe Cr and 
Ni-Co-Cr Mo alloys. Data, where available, were taken from the literature. 
However, values for y c and D c c were determined by fitting measured and 
predicted C profiles in binary alloys which were carburized in an atmosphere 
which did not form carbides. Thermodynamic data for the more complex 
carbides were determined by systematically fitting measured and predicted 
C profiles in the progressively more complex alloys. The wide range of 
solubility of some carbides (i.e., C^Q and Cr 7 C. 3 , can dissolve significant 
amounts of Fe, Ni and Mo) was approximated by several carbides with 
selected discrete stoichiometries. For example, in Ni Fe Cr alloys, the Cr/Fe 
ratio changes continuously in the M7C3 carbide. This carbide was therefore 
approximated by the two carbides of fixed stoichiometry, Cr 4 Fe-?C 3 and 
Cr 5 Fe 2 Ci. The equilibrium constants were determined by assuming the 
mixed Fe Cr carbides to be regular solutions of the stable phase Cr 7 C}. The 
model gave very good agreement with the predicted C concentration profiles 
for the binary Ni Cr alloys and also agreed well with the new experimental 
results for the carbide distributions. In the ternary and quaternary alloys, 
reasonable agreement was also found considering the simplifying assumption 
of approximating the continuously varying carbide compositions with 
discrete stoichiometric carbides. 

Farkas and Ohla 14 followed the two-step approach of Bongartz et a/. 18 
in modeling the simultaneous precipitation of up to three carbides but 
included the ternary cross-term effect of a substitutional element. A near- 
surface region was depleted of Cr in three commercial superalloys and two 
Fe Ni Cr alloys by preoxidizing at 900 C for 100 h. Similar to the studies 
of Bongartz, measured and predicted C profiles were used to derive approxi- 
mate values for D cx in the alloys. The derived values for D cx were in 
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Fig. 1 1 . Measured and predicted C concentration pro- 
files for the carburization of a Ni 30Cr alloy. (After 
Engstrom et at.' 2 ) 


reasonable agreement with values determined in other studies. The measured 
and predicted C profiles showed the need to account for the ternary cross- 
term effect when predicting carburization in multicomponent alloys. Farkas 
et aI M later applied a similar model to the carburization and sulfidation of 
Ni Cr Mo alloys. As expected, for the same time and temperature condi- 
tions, S penetration and sulfide precipitation are significantly less than that 
associated with carburization since S diffuses as a substitutional element in 
contrast to the interstitial diffusion of C. 

Engstrom et al. y2 also followed the two-step approach of Bongartz but 
developed a generalized numerical model which makes use of a thermo- 
dynamic package to make all equilibrium calculations. There are no restric- 
tions on the number of components or phases which occur in the problem 
to be modeled as long as thermodynamic and kinetic data are available. The 
model also utilizes a labyrinth factor equal to the square of the volume 
fraction of remaining matrix. This model, which incorporates no adjustable 
parameters, was also applied to the carburization of Ni Cr alloys and eval- 
uated with previously published experimental data. 24 '™ Because of the direct 
link with the thermodynamic program, the model was able to account for 
the variable composition of the carbides which were modeled discretely by 
Bongartz. 24 All other data used by the model w ; ere determined independently. 
There was reasonable agreement between the calculated C profiles and car- 
bide distributions, as shown in Fig. 11. Although only data measured 
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independently were used in the calculations, the values used (e.g., D c x ) were 
not reported to allow a comparison to the values determined in the progress- 
ive manner by Bongartz. 29 

The second method to model the carburization process when accom- 
panied by carbide precipitation was presented by Sockel cl air 8 In this 
approach, the precipitation reaction is included in FicVs second law. Follow- 
ing the treatment by Sockel, the total C concentration is given as 
Cc >t = C(^ + Cc, where C ppt is the precipitated C in the carbide and C c is 
the C in solution. The concentration term in Fick’s second law, PC/ f t (Eq. 
(2a)), must be modified to include the total C concentration, C[ ot , w r hich 
will change with time. However, other terms in Fick's second law remain 
unchanged since the precipitated C does not contribute to the C flux. By 
assuming ideal solution behavior, the solubility product for a B b C a carbide 
is given as £ v =CbCc , the concentration of precipitated C becomes 


C7 pt = 


a 

h 



Kl 

Cc 


(17) 


where C is the total concentration of B at a particular location in the 
carbide subscale layer, defined as C£ pl + G*. Inserting these equations into 
Fick's second law, and assuming that D c , c is independent of concentration, 
yields 


1 + 


K s 


b) \Ci 


ft PX 2 


C c 


(18) 


F-D equivalents can again be substituted for the derivatives in this equation, 
however, the equation is now ; nonlinear and the concentration of C in solu- 
tion, Cc, at the new' time f + A/ cannot be solved for directly but requires 
an iterative numerical technique. 9 Sockel 28 utilized these equations for the 
case where the solute B was immobile w f hich yields a constant value for 
C B >1 in Eq. (17). Hence, when C c was determined in Eq. (18), the value for 
C ppt could be calculated using Eq. (17) and C B and the totals, C{° 1 and 
C B 1 could be easily determined. If B is mobile, two equations similar to Eqs. 
(17) and (18) could be derived also for B. However, for each node within 
the carbide subscale zone, the F-D equivalents to Eq. (18) for both B and 
C would have to be each iteratively solved to obtain values for C B and C\ 
which also satisfied the solubility product, a significant amount of calculation 
for the concentrations at each node! 

Sockel el a/. 28 utilized the approach outlined above to model the carburi- 
zation of Ni 25Cr alloys under conditions where Cr^C?, Cr 7 C\ and Cr 27 C b 
form. A labyrinth factor, approximately equal to the volume fraction of the 
remaining matrix, was used. Sockel assumed a constant value for Dc , an 
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ideal solution model, and an immobile Cr solute. Data were taken from the 
literature for D t and K s for each of these carbides, these latter values which 
had to be “corrected” in order to be used in the computer calculations. 
Good agreement was obtained between the measured and predicted weight 
gain during carburization of Ni 25Cr alloys. There was also reasonable 
agreement between the measured and predicted volume fraction of carbides 
versus distance in these alloys. As expected, the predictions also show that 
the thickness of each of the carbide layers increases parabolically with time. 


FUTURE 

The intent of the preceding discussion was to demonstrate that F-D 
techniques have the flexibility to simulate the diflfusional transport associated 
with most high-temperature corrosion problems of interest. The tremendous 
advances in computing power have taken the solution of these problems 
from the exclusive domain of mainframe computers to the common task of 
desktop PCs in less than 10 years. However, as at the time of the earliest 
studies using numerical methods to model diffusional transport, the chief 
impediment in performing numerical modeling of diffusion processes lies in 
the lack of accurate thermodynamic (e.g., phase stabilities and solubilities) 
and kinetic data (e.g., diffusion coefficients). Most of the studies discussed 
above cite this as an obstacle to more detailed analyses. Unfortunately, one 
of the limitations of the F-D technique is that it cannot be “inverted” to 
directly yield concentration-dependent diffusivity data, although several of 
the studies discussed above used the technique to “back-calculate” apparent 
diffusivities (see Refs. 18, 29, 30). 

Two statements can be made here for future directions. First, is the 
need for more work in line with the approach taken by Engstrom 32 in linking 
predictive thermodynamic packages with diffusion models. With sufficient 
data, these thermodynamic programs can predict the phase relationships 
and boundaries for multiphase and multicomponent alloys critical to the 
application of the diffusion models. However, a dearth of diffusivity data, 
especially in multicomponent alloys, continues. This lack of data is not likely 
to be quickly remedied considering the extensive amount of concentration 
data required to measure Ds in ternary and higher order systems. Morral, 
Thompson and co-workers, 33 ’ 34 have developed a different approach of pre- 
dicting diffusion in multicomponent, multiphase alloys. This approach util- 
izes a “square-root diffusivity” matrix [r] to predict interdiffusion. Although 
the analysis assumes concentration independence, the elements of the [r] 
matrix can be determined with less experimental work than the traditional 
Boltzmann Matano analysis to yield [/)]. In any multicomponent system, 
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the [r] matrix is related to the familiar diflusivities by [£>] = [?*][/]. The 
accuracy and usefulness of the approach for multicomponent alloys has been 
demonstrated in predicting concentration profiles in Ni-Cr-Al-Mo alloys . 35 
The approach has also been applied to predict coating-substrate interdiffu- 
sion in multicomponent systems ' 36 and to investigate the effect of third ele- 
ment additions in helping to establish an external, protective, solute scale . 37 
Incorporating the square root diffusivity analysis into the F-D technique 
may allow modeling diffusional transport associated with high-temperature 
corrosion in more complex commercial alloys. Certainly the necessary kinetic 
data, nearly impossible to measure using traditional approaches, will be 
easier to measure with this approach . 38 

Lastly, most diffusion analyses are performed on relatively simple geo- 
metries, planar, cylindrical or spherical. In applications where the shape 
cannot be reduced to these simple geometries, more complex, commercial, 
heat transfer codes may be used. In this case, the temperature-dependent 
thermal diffusivity is replaced by a concentration-dependent kinetic diffusiv- 
ity. These codes, such as SINDA 39 (Systems Improved Numerical Differenc- 
ing Analyzer), typically have considerable flexibility in handling time- 
dependent boundary conditions, complex starting conditions, and moving 
boundaries. The drawback to these codes is that they are often extensive 
programs containing numerous capabilities and are generally resident on 
mainframe computers. 
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